14. Análisis de Asociación de Variables Aleatorias

En esta sección, se presentan tests para analizar

• la asociación entre una v.a. continua y una v.a. categórica (ANOVA), o

• la asociación entre variables aleatorias discretas o categóricas (Chi-cuadrado, Bowker).

14.1. Repaso (prueba de hipótesis)

1) Revisamos la prueba de entrada sobre “Inferencia en Estadística Paramétrica”

2) ¿Cómo reportamos la prueba de t?

Se realizó una prueba t [de dos muestras o pareada] para comparar [variable de respuesta de interés] en [grupo 1] y [grupo 2]. Hubo [o no hubo] una diferencia significativa en [variable de respuesta de interés] entre [grupo 1] (M = [media], SD = [desviación estándar]) y [grupo 2] (M = [media], sd = [desviación estándar]), t(df) = [valor t], p = [valor p].

Pero hay flexibilidad en reportar la prueba. Por ejemplo:

Se realizó una prueba t [de dos muestras o pareada] para comparar [variable de respuesta de interés] en [grupo 1] y [grupo 2]. El [grupo 1] tuvo [variable de respuesta de interés] mayor [o menor] (M = [media], SD = [desviación estándar]) que los del [grupo 2] (M = [media], SD = [desviación estándar]), t(df) = [valor t], p = [valor p].

3) ¿Cómo es la tabla de decisión vs. estado real? ¿Qué es error de tipo I, nivel de significancia del test, potencia?

4) ¿Cómo es el test de medias para dos poblaciones normales con la misma varianza desconocida?

Test de medias para dos poblaciones normales con la misma varianza desconocida (caso común)

Sean \(X_1,\cdots,X_n\) e \(Y_1,\cdots,Y_m\) muestras independientes de poblaciones normales con medias desconocidas \(\mu_x\) y \(\mu_y\) y misma varianza desconocida \(\sigma^2\). Consideremos el test de hipótesis:

\[H_0: \mu_x = \mu_y\]

\[H_1: \mu_x \neq \mu_y\]

del Corolario del Teo de Fisher-Cochran se cumple:

\[ \frac{\bar{X} - \bar{Y} - (\mu_x - \mu_y)}{\sqrt{S_p^2(\frac{1}{n}+ \frac{1}{m})}} \sim t_{n+m-2}\]

donde

\[S_p^2 = \frac{(n-1)S_X^2 + (m-1)S_Y^2}{n+m-2}\]

de manera que se rechaza \(H_0\) si

\[\frac{\bar{X} - \bar{Y}}{\sqrt{S_p^2(\frac{1}{n}+ \frac{1}{m})}} > t_{\frac{\alpha}{2},n+m-2}\]

14.2. ¿Por qué ANOVA?

¿Qué occure si tenemos 3, 5, o 10 grupos/modelos/métodos para comparar? ¿Qué problemas tendremos al utlizar t-tests?

1. En el caso en que se requiera comparar más de 2 grupos, el t-test se vuelve ineficiente, porque no queremos hacer un montón de t-tests para cada par.

2. Además, family-wise error rate (la tasa de error de la familia, error global) que es la probabilidad de cometer al menos un error de Tipo I en múltiples pruebas estadísticas realizadas en los mismos datos aumenta. Para \(c\) tests se calcula como:

\[\alpha_{global} = 1-(1-\alpha)^c\]

• con c=2 tests, el tasa de error de tipo I es 0.0975 (alrededor de 2*0.05=0.1)

• con c=3 tests, el tasa de error de tipo I es 0.143 (alrededor de 3*0.05=0.15)

• con c=10 tests, el tasa de error de tipo I es 0.40 (alrededor de 10*0.05=0.5)

Por lo tanto, es mejor abordar con el modelo del Análisis de la Varianza (ANOVA), o también llamado análisis de factores, para estudiar el efecto de uno o más factores (cada uno con dos o más niveles), o para comparar más de 2 grupos, sobre la media de una variable continua.

Si queremos hacer la comparación al mismo tiempo, ¿cómo sería el H0 si tenemos \(m\) grupos?

\[H_0: \mu_1 = \mu_2 =\cdots =\mu_m\qquad\]

\[H_1: \exists i,j \in \{1,\cdots m \},\, \mu_i \neq \mu_j\]

¿Qué distribución podemos utilizar?

14.3. Repaso (chi-cuadrado)

14.3.1. La distribución chi-cuadrado

Sean \(Z_1,\cdots, Z_n\, v.a. i.i.d. \, \sim {\it N}(0,1)\) entonces

\[ Y = Z_1^2+\cdots+Z_n^2 \sim \chi_{n}^2\]

donde \(n\) son los grados de libertad de la distribución.

Propiedades de la distribución \(\chi^2\):

(i) Propiedad aditiva:

Si \(X_1\) y \(X_2\) son dos v.a. independientes distribuidas \(\chi^2\) de \(n_1\) y \(n_2\) grados de libertad, entonces

\(X_1+X_2 \sim \chi_{n_1+n_2}^2\) (importante)

(ii) Esperanza y Varianza:

\( E[X]= n\) (importante)

\(Var[X]= 2n\)

suppressMessages(library(dplyr))
suppressMessages(library(plotly))
suppressMessages(library(ggplot2))
suppressMessages(library(rmarkdown))
vec <- seq(0,100,by=0.1)
params <- c(1:5, seq(10, 100, by=10)) 
pvec <- list()
for (i in 1:length(params)){
    pvec[[i]] <- dchisq(vec,df=params[i],ncp=0)
}
steps <- list()
fig <- plot_ly(width=600,height=600) %>% layout(title = "\n \n Densidad de Probabilidad Chi-cuadrado", yaxis = list(range=c(0,0.3)))
for (i in 1:length(params)){
    fig <- add_lines(fig, x=vec, y=pvec[[i]], visible=if (i==1) TRUE else FALSE, mode='lines', line=list(color='blue'), showlegend=FALSE)
    steps[[i]] = list(args = list('visible', rep(FALSE, length(params))), label=params[i], method='restyle')
    steps[[i]]$args[[2]][i] = TRUE
}
fig <- fig %>% layout(sliders = list(list(active=0, currentvalue = list(prefix = "df: "), steps=steps)))
fig

¿Qué sigue la distribución chi-cuadrado?

14.3.2. La distribución de la varianza (y media) muestral del caso Normal

Teorema de Fisher-Cochran

Sean \(X_1,\cdots,X_n\) v.a. i.i.d. \({\cal N}(\mu,\sigma^2)\) entonces la media y varianza muestral cumplen:

\(\begin{equation} \begin{array}{lcll} (i) & \bar{X} &\sim& {\cal N}(\mu, \frac{\sigma^2}{n})\\ \\ (ii) & {\displaystyle \frac{(n-1)S^2}{\sigma^2}}& \sim& \chi_{n-1}^2 \text{ (<== enfoque en esta clase)}\\ \\ (iii)& \bar{X} &{\mathrel \perp} & S^2 \quad \text{(independentes)}\\ \end{array} \end{equation}\)

\(S^2\) es la varianza muestral. ¿Cuál es la diferencia con \(\sigma^2\)? Cómo calculamos \(S^2\)?

\[S^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i-\bar{X})^2 \]

Por qué en \( {\displaystyle \frac{(n-1)S^2}{\sigma^2}} \sim \chi_{n-1}^2\), el grado de liberdad es \(n-1\)?

¿Cómo entendemos la segunda propiedad intuitivamente?

\(\begin{equation} \begin{array}{rll} {\displaystyle \frac{(n-1)\frac{1}{n-1}\sum_{i=1}^n (X_i-\bar{X})^2}{\sigma^2}} & \sim& \chi_{n-1}^2 \\ {\displaystyle \frac{\sum_{i=1}^n (X_i-\bar{X})^2}{\sigma^2}} & \sim & \chi_{n-1}^2 \\ \sum_{i=1}^n ({\displaystyle \frac{X_i-\bar{X}}{\sigma}} )^2 & \sim & \chi_{n-1}^2 \end{array} \end{equation}\)

¿Qué pasa si dividimos un chi-cuadrado por otro chi-cuadrado?

14.4. La distribución F

Sean \(X \sim \chi_n^2\) e \(Y \sim \chi_m^2\) dos v.a. independientes \(\chi^2\) de grados de libertad \(n\) y \(m\) respectivamente, entoncese se define:

\[F = \frac{\frac{X}{n}}{\frac{Y}{m}} \sim F_{n,m}\]

donde \(F_{n,m}\) es la distribución \(F\) de \(n\) y \(m\) grados de libertad. También se nota \(F(n,m)\).

vec <- seq(0,5,by=0.01)
params <- seq(1,20,by=1)
pvec <- list()
for (i in 1:length(params))
    for (j in 1:length(params)){
        k = length(params)*(i-1) + j
        pvec[[k]] <- df(vec,df1=params[i],df2=params[j],ncp=0)
}

steps1 <- list()
steps2 <- list()
fig <- plot_ly(width=600,height=600) %>% layout(title = "\n \n Densidad de Probabilidad F(n, m)",
                                                 yaxis = list(range=c(0,1)))
for (i in 1:length(params)){
     for (j in 1:length(params)){
        k = length(params)*(i-1) + j
        fig <- add_lines(fig, x=vec, y=pvec[[k]], 
                     visible=if ((i==1) && (j==1)) TRUE else FALSE,
                     mode='lines', line=list(color='blue'), showlegend=FALSE)
        steps2[[j]] = list(args = list('visible', rep(FALSE, length(params)*length(params))), 
                      label=params[j], method='restyle')
        steps2[[j]]$args[[2]][k] = TRUE
        steps1[[i]] = list(args = list('visible', rep(FALSE, length(params)*length(params))), 
                      label=params[i], method='restyle')
        steps1[[i]]$args[[2]][k] = TRUE
  }                   
}
fig <- fig %>% layout(sliders = 
                      list( list(active=0, currentvalue = list(prefix = "df1 (n): "), pad = list(t=20), steps=steps1),
                      list(active=0, currentvalue = list(prefix = "df2 (m): "), pad = list(t=100), steps=steps2)))
fig

Percentil

Sea \(F_{\alpha,n,m}\) tal que \(P\{F > F_{\alpha,n,m} \} = \alpha\), \(F_{\alpha,n,m}\) es el percentil \(100(1-\alpha)\) de la distribución \(F_{n,m}\).

Relación entre las distribuciones F y t

Si \(X \sim t_{n}\) entonces \(X^2 \sim F_{1,n}\). Por qué?

Por construcción, si \(X \sim t_{n}\), entonces

\[X = \frac{Z}{\sqrt{\frac{Y}{n}}}, \qquad con \, Z \sim {\cal N}(0,1) \,e \,Y \sim \chi_n^2\]

Entonces

\[X^2 = \frac{\frac{Z^2}{1}}{\frac{Y}{n}}, \qquad con \, Z^2 \sim \chi_1^2 \,e \,Y \sim \chi_n^2\]

14.5. ANOVA de un factor (unidireccional, one-way)

Tenemos \(H_0\) y \(H_1\) en las palabras así:

\[H_0: \text{Todas las medias de los grupos son iguales / no hay diferencias en las medias de los grupos}\]

\[H_1: \text{Al menos dos medias de los grupos son diferentes / las medias de los grupos no son iguales}\]

Formalmente, sean \(m\) muestras (grupos) independientes de tamaño \(n\):

\[X_{i1},\cdots,X_{in} \sim i.i.d. \,{\cal N}(\mu_i,\sigma^2), \qquad i=1,\cdots, m\]

donde \(\mu_i\) y \(\sigma^2\) son desconocidos. Estamos interesados en probar:

\[H_0: \mu_1 = \mu_2 =\cdots =\mu_m\qquad\]

\[H_1: \exists i,j \in \{1,\cdots m \},\, \mu_i \neq \mu_j\]

Ahora veamos un contexto concreto. Suponga que tenemos \(m\) tratamientos distintos, donde el resultado de aplicar el tratamiento \(i\) a un individuo es una v.a. normal. Estamos interesados en probar la hipótesis de que todos los tratamientos tienen el mismo efecto. Para ello se aplica cada tratamiento a una muestra distinta de \(n\) individuos y se analizan los resultados.

Decimos que ANOVA es para analizar la asociación entre una v.a. continua (variable dependiente) y una v.a. categórica (variable independiente). En este ejemplo, ¿cuál es la v.a. continua? ¿Cuál es la v.a. categórica? Cuál es el factor? ¿Puedes formular las hipótesis de otra manera en este ejemplo como “… (no) depende de …”?

\(H_0\): El resultado de aplicar el tratamiento no depende de (las categorías del) tratamientos.

\(H_1\): El resultado de aplicar el tratamiento depende de (las categorías del) tratamientos.

14.5.1. Supuestos (y robustez)

¿Qué supuestos tiene el ANOVA según las definiciones matemáticas anteriores?

Consideramos sus supuestos y robustez aquí:

• Normalidad: cada muestra/grupo (o v.a. \(X_{ij}\)) debe provenir de una distribución normal. Utilizar por ejemplo, test de Shapiro-Wilk. Utilizar Q-Q plot para descartar la presencia de outliers.

  • El ANOVA unidireccional se considera una prueba robusta contra el supuesto de normalidad, cuando el tamaño de cada muestra (o grupo) no es pequeño (y los tamaños de grupos son iguales); si no, consideremos test nonparamétrico como Kruskal-Wallis H Test.

• Varianza común (homogeneidad): misma varianza en los grupos. Se puede probar la igualdad de varianzas con el test de Levene o con el test de Brown-Forsythe.

  • Si este supuesto falla, consideraremos test nonparamétrico como Kruskal-Wallis H Test.

• Observaciones independientes. Las poblaciones de los grupos son independientes y las observaciones de cada grupo son independientes. Este supuesto no puede ser probado con ningún estadístico, es una consideración de diseño

  • La violación de esto es muy grave (por ejemplo, la tasa de error de Tipo I está sustancialmente inflada).

14.5.2. Idea general de la derivación ANOVA

Hay dos puntos fundamentos para entender la derivación ANOVA:

1) varianza intragrupos (within-group) vs. varianza entregrupos (between-group)

(fuente)

2) Formalmente, construimos dos estimadores de la varianza común \(\sigma^2\), y comparar estas dos “varianzas”

• El primer estimador no depende de si \(H_0\) es cierto o no.

• En cambio, el segundo estimador asume que \(H_0\) es cierto (en caso contrario este estimador sobreestima a \(\sigma^2\)).

El test compara ambos estimadores por una división entre ellos: \({\displaystyle \frac{\text{segundo estimador}}{\text{primer estimador}}}\).

Rechazamos \(H_0\) cuando la tasa entre el segundo y primer estimador es suficientemente grande.

Si los dos estimadores son sobre varianza intragrupos y varianza entregrupos, cómo construimos ellos? Elige una opción:

A. \({\displaystyle \frac{\text{segundo estimador}}{\text{primer estimador}} = \frac{\text{varianza intragrupos}}{\text{varianza entregrupos}}}\)

B. \({\displaystyle \frac{\text{segundo estimador}}{\text{primer estimador}} = \frac{\text{varianza entregrupos}}{\text{varianza intragrupos}}}\)

14.5.3. El primer estimador de \(\sigma^2\) que no depende de \(H_0\) (\(SS_W\))

Definimos:

\[X_{i.} = \frac{\sum_{j=1}^n X_{ij}}{n}\]

Recordemos la propiedad:

\[{\displaystyle \frac{(n-1)S^2}{\sigma^2}} \sim \chi_{n-1}^2\]

Así que:

\[ \sum_{j=1}^n \frac{(X_{ij} - X_{i.})^2}{\sigma^2} \sim \chi^2_{n-1} \quad \text{ (Por qué?)}\]

\[\sum_{i=1}^m \sum_{j=1}^n \frac{(X_{ij} - X_{i.})^2}{\sigma^2} \sim \chi^2_{nm-m} \quad \text{ (Por qué?)}\]

Denominaremos

\[ SS_W = \sum_{i=1}^m \sum_{j=1}^n (X_{ij} - X_{i.})^2\]

(\(SS_W\) es “within samples sum of squares”, o “residual sum of squares”)

Sustituimos la parte izquierda de la distribución chi-cuadrado anterior:

\[ \frac{SS_W}{\sigma^2} \sim \chi^2_{nm-m}\]

Entonces

\[ \frac{E[SS_W]}{\sigma^2} = nm-m \quad \text{ (Por qué?)}\]

y luego

\[ \frac{E[SS_W]}{nm-m} = \sigma^2\]

con lo cual tenemos nuestro primer estimador insesgado de \(\sigma^2\):

\[MSE = \frac{SS_W}{nm-m}\]

¿Por qué insesgado?

Nota

Para un parámetro \(\theta\) y un estimador \(\hat \theta\), si

\[E[\hat \theta] = \theta\]

entonces \(\hat \theta\) es un estimador insesgado de \(\theta\). Es decir, la distribución de un estimador insesgado se centra en el parámetro verdadero.

¿\({\displaystyle MSE = \frac{SS_W}{nm-m}}\) corresponde a la varianza entregrupos o la varianza intragrupos?

\({\displaystyle MSE = \frac{SS_W}{nm-m}}\) se conoce como “error cuadrático medio” (mean square error), “cuadrado de la media de los residuos” (residual mean square) o “varianza intragrupos” (within-group variance). Notar que no considera ningun supuesto sobre \(H_0\).

14.5.4. El segundo estimador de \(\sigma^2\) que depende de \(H_0\) (\(SS_B\))

En este caso suponemos que se cumple \(H_0\), es decir

\[\mu = \mu_1 = \mu_2 = \cdots = \mu_m\]

De manera de que las medias muestrales cumplen:

\[X_{1.},X_{2.},\cdots,X_{m.} \sim i.i.d. \, {\cal N}\left(\mu, \frac{\sigma^2}{n} \right) \quad \text{ (Por qué?)}\]

entonces:

\[ \frac{\sum_{i=1}^m (X_{i.}-X_{..})^2}{\frac{\sigma^2}{n}} \sim \chi^2_{m-1} \quad \text{ (Por qué?)}\]

Cómo aplicamos la segunda propiedad del Teorema de Fisher-Cochran aquí?

Nota

Sean \(X_1,\cdots,X_n\) v.a. i.i.d. \({\cal N}(\mu,\sigma^2)\) entonces la varianza muestral cumplen:

\(\begin{equation} \begin{array}{lcll} (ii) & {\displaystyle \frac{(n-1)S^2}{\sigma^2}}& \sim& \chi_{n-1}^2\\ \iff & {\displaystyle \frac{(n-1)\frac{1}{n-1}\sum_{i=1}^n (X_i-\bar{X})^2}{\sigma^2}} & \sim& \chi_{n-1}^2 \\ \iff & {\displaystyle \frac{\sum_{i=1}^n (X_i-\bar{X})^2}{\sigma^2}} & \sim & \chi_{n-1}^2 \end{array} \end{equation}\)

Denominaremos

\[ SS_B = n\sum_{i=1}^m (X_{i.} - X_{..})^2 = \sum_{i=1}^m \sum_{i=1}^n (X_{i.} - X_{..})^2\]

(\(SS_B\) es “between samples sum of squares”, o “model sum of squares”.)

Sustituimos la parte izquierda de la distribución chi-cuadrado anterior:

\[ \frac{SS_B}{\sigma^2} \sim \chi^2_{m-1}\]

Entonces, si \(H_0\) es cierto, se cumple

\[ \frac{E[SS_B]}{\sigma^2} = m-1 \quad \text{ (Por qué?)}\]

y luego

\[ \frac{E[SS_B]}{m-1} = \sigma^2\]

con lo cual tenemos nuestro segundo estimador insesgado de \(\sigma^2\):

\[MSB = \frac{SS_B}{m-1}\]

(¿Por qué insesgado?)

¿\({\displaystyle MSB = \frac{SS_B}{m-1}}\) corresponde a la varianza entregrupos o la varianza intragrupos?

\({\displaystyle MSB = \frac{SS_B}{m-1}}\) se conoce como “cuadrado de la media entre muestras” (between samples mean square), “cuadrado de la media del modelo” (model mean square) o “varianza entregrupos” (between-group variance).

Ahora, tenemos todos para definie un Test estadístico!

14.5.5. Prueba F (F-test)

Se define como el estadístico de prueba (Test Statistic):

\[TS = \frac{\text{MSB}}{\text{MSE}} = \frac{\text{varianza entregrupos}}{\text{varianza intragrupos}}\]

Cuándo rechazamos la hipótesis nula?

Se rechazará \(H_0\) cuando TS sea suficientemente grande.

¿Cuál es la distribución de TS?

\[TS = \frac{MSB}{MSE} = \frac{\displaystyle \frac{SS_B}{m-1}}{\displaystyle \frac{SS_W}{nm-m}}\]

Sabemos:

(1) Dos distribuciones chi-cuadrado:

\[ \frac{SS_B}{\sigma^2} \sim \chi^2_{m-1}\]

\[ \frac{SS_W}{\sigma^2} \sim \chi^2_{nm-m}\]

(2) La definicizón de la distribución F:

Sean \(X \sim \chi_n^2\) e \(Y \sim \chi_m^2\) dos v.a. independientes \(\chi^2\) de grados de libertad \(n\) y \(m\) respectivamente, entoncese se define:

\[F = \frac{\frac{X}{n}}{\frac{Y}{m}} \sim F_{n,m}\]

Así que:

\[ TS \sim F_{m-1,nm-m}\]

Utilizaremo F en lugar de TS para referir al estadístico de prueba desde ahora:

\[ F \sim F_{m-1,nm-m}\]

¿Cómo utilizamos p-value o valor crítico para hacer la prueba de hipótesis? (Se rechazará \(H_0\) cuando F sea suficientemente grande.)

1) p-value

Se calcula \(f\) en una muestra dada, y

\[p = P_{H_0}\{F > f\} \]

¿Cuándo rechazamos \(H_0\)?

Si \(p<\alpha\), se rechaza \(H_0\). No se rechaza \(H_0\) en caso contrario.

2) Valor crítico

Sea entonces \(F \sim F_{m-1,nm-m}\) y \(F_{\alpha,m-1,nm-m}\) tal que:

\[P_{H_0}\{F > F_{\alpha,m-1,nm-m}\} = \alpha\]

Se calcula \(f\) en una muestra dada.

¿Cuándo rechazamos \(H_0\)?

Si \(f > F_{\alpha,m-1,nm-m}\), se rechaza \(H_0\). No se rechaza \(H_0\) en caso contrario.

14.5.6. Identidad de suma de cuadrados

¿Cómo calculamos la varianza muestral de la muestra \(\{X_{ij}\}\) (i=1, …m; j=1,…n)?

Definimos \(SS_{total} = SS_{T}\) (total sum of squares) como:

\[\sum_{i=1}^m \sum_{j=1}^n (X_{ij} - X_{..})^2\]

Cumple:

\[\sum_{i=1}^m \sum_{j=1}^n (X_{ij} - X_{..})^2 = \sum_{i=1}^m \sum_{j=1}^n (X_{ij} - X_{i.})^2 + \sum_{i=1}^m \sum_{j=1}^n (X_{i.} - X_{..})^2 \quad \text{ (Por qué?)}\]

Así que tenemos:

\[SS_{total} = SS_{intragrupos(within)} + SS_{entregrupos(between)}\]

\[SS_{T} = SS_{W} + SS_{B}\]

14.5.7. Ejemplo (industria de automóviles)

Una industria de automóviles desea probar la calidad de 3 tipos de combustibles. Para ello observa el rendimiento (millas por 10 galones de combustible) de los 3 tipos de combustibles en 15 vehículos idénticos (5 con cada tipo de combustible). (Asumimos que los supuestos de ANOVA se mantienen aquí, pero hay que examinar los supuestos en tu análisis de los datos.)

Cuál es la hipótesis nula aquí?

Veamos primero la distribución de los datos:

library(tidyr)
library(ggplot2)
library(dplyr)
library(cowplot)

options(repr.plot.width=10, repr.plot.height=5)
options(digits=4)

gas1 <- c(220,251,226,246,260)
gas2 <- c(244,235,232,242,225)
gas3 <- c(252,272,250,238,256)

df <- do.call(rbind, Map(data.frame, gas1=gas1, gas2=gas2, gas3=gas3))
df <- df %>% pivot_longer(cols=c('gas1', 'gas2', 'gas3'),names_to='combustible', values_to='milla')
df

df %>% group_by(combustible) %>% summarise(M=mean(milla), SD=sd(milla))

#mean_sdl: computes the mean plus or minus a constant times the standard deviation. By default mult=2
p1 <- ggplot(df, aes(x=combustible, y=milla)) + 
    geom_violin(trim=FALSE) + 
    geom_dotplot(binaxis='y', stackdir='center', dotsize=1, binwidth=1.5) +
    stat_summary(fun.data=mean_sdl, fun.args=list(mult=1), geom="pointrange", color="red", lwd=1) + 
    theme(text=element_text(size=14))
p2 <- ggplot(df, aes(x=combustible, y=milla)) + geom_boxplot() + theme(text=element_text(size=14))

plot_grid(p1, p2, align="h", nrow=1, ncol=2, rel_widths= c(1/2, 1/2))

A tibble: 15 x 2

combustible	milla
<chr>	<dbl>
gas1	220
gas2	244
gas3	252
gas1	251
gas2	235
gas3	272
gas1	226
gas2	232
gas3	250
gas1	246
gas2	242
gas3	238
gas1	260
gas2	225
gas3	256

A tibble: 3 x 3

combustible	M	SD
<chr>	<dbl>	<dbl>
gas1	240.6	16.965
gas2	235.6	7.701
gas3	253.6	12.280

En tu tarea, tienes que desarrollar el test. Aquí, mostramos el resultado de R:

# Método 1
m <- aov(milla ~ combustible, data=df)
summary(m)

            Df Sum Sq Mean Sq F value Pr(>F)
combustible  2    863     432     2.6   0.12
Residuals   12   1992     166               

Otra manera de hacer ANOVA en R:

# Método 2
m0 <- lm(milla ~ combustible, data=df)
anova(m0)

A anova: 2 x 5

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
	<int>	<dbl>	<dbl>	<dbl>	<dbl>
combustible	2	863.3	431.7	2.601	0.1152
Residuals	12	1991.6	166.0	NA	NA

14.5.8. Caso de muestras con diferentes tamaños

¿Qué ocurre si las muestras no son del mismo tamaño?

\[\sum_{i=1}^m \sum_{j=1}^{n_i} \frac{(X_{ij} - X_{i.})^2}{\sigma^2} \sim \chi^2_{N} \quad \text{ (Por qué?)}\]

donde \(N = \sum_{i=1}^m n_i - m\). Y entonces,

\[ SS_W = \sum_{i=1}^m \sum_{j=1}^{n_i} (X_{ij} - X_{i.})^2\]

por otra parte

\[\sum_{i=1}^m \frac{(X_{i.}-X_{..})^2}{\frac{\sigma^2}{n_i}} \sim \chi^2_{m-1}\]

y se definen

\[ SS_B = \sum_{i=1}^m n_i(X_{i.}-X_{..})^2 = \sum_{i=1}^m \sum_{j=1}^{n_i} (X_{i.}-X_{..})^2\]

\[ F = \frac{\frac{SS_B}{m-1}}{\frac{SS_W}{N}} \sim F_{m-1,N} \qquad (N= \sum_{i=1}^m n_i - m)\]

Finalmente, si \(F > F_{\alpha,m-1,N}\) se rechaza \(H_0\). O también puedes utilizar p-value.

14.6. Métodos de comparaciones múltiples

Si el test ANOVA es estadísticamente significativo, lo único que podemos concluir es que una o mas de las diferencias entre grupos es significativa, pero no sabemos que grupos son los que difieren.

Para ello es necesario realizar “comparaciones múltiples”. Se trata de comparar cada grupo con otro grupo, o un promedio de grupos se puede comparar con otros. Se consideran dos posibilidades:

• Comparaciones planeadas: existe interés por algunos grupos en particular

• Comparaciones post hoc: no existen hipótesis específicas de algunos grupos. Nos centramos en este caso aquí.

Aquí presentamos algunos métodos de comparaciones múltiples que son los más comúnmente usados.

¿Recuerdas cómo calculamos la tasa de error de tipo I si hacemos \(c\) comparación por pares en las que cada uno utlizamos \(\alpha=0.05\) en cada prueba de hipótesis?

Nota

La tasa de error familiar / error global (family-wise error rate) es la probabilidad de cometer al menos un error de tipo I en múltiples pruebas estadísticas realizadas en los mismos datos. Para \(c\) tests se calcula como:

\[\alpha_{global} = 1-(1-\alpha)^c\]

Si utilizamos \(\alpha=0.05\) en cada prueba de hipótesis:

• con c=2 tests, el tasa de error global es 0.0975, alrededor de 2*0.05=0.1

• con c=10 tests, el tasa de error global es 0.40, alrededor de 10*0.05=0.5

Es decir, con c=10 tests, \(\alpha=0.05\), la probabilidad de cometer al menos un error de tipo I (\(\alpha_{global}\)) es alrededor de 0.5.

Con c tests, la probabilidad de cometer al menos un error de tipo I (\(\alpha_{global}\)) es alrededor de \(c\alpha\).

¿Cómo podemos modificar \(\alpha\) (de cada prueba) para que \(\alpha_{global}\) se mantenga 0.05?

14.6.1. Corrección de Bonferroni

Se calcula un nuevo valor de \(\alpha\) para mantener el error de tipo I global (familiar):

\[ \alpha_{nuevo} = \frac{\alpha}{c}\]

Por ejemplo, si \(c=10\) y \(\alpha=0.05\), entonces \(\alpha_{nuevo} = \frac{0.05}{10} = 0.005\)

Ejemplo (fuente)

Supongamos que un profesor quiere saber si tres técnicas de estudio diferentes conducen a puntuaciones diferentes en los exámenes. Para comprobarlo, asigna aleatoriamente a 30 estudiantes cada técnica de estudio. Después de una semana de utilizar la técnica de estudio asignada, cada estudiante realiza el mismo examen.

(Asumimos que los supuestos de ANOVA se mantienen aquí, pero hay que examinar los supuestos en tu análisis de los datos.)

1) Realiza un ANOVA unidireccional. ¿Cuál es \(H_0\)?

2) Descubre que el valor p (de prueba F de ANOVA) es 0.0476. ¿Se rechaza \(H_0\) o no?

Como p<0.05, rechaza \(H_0\) y concluye que no todas las técnicas de estudio producen la misma puntuación media en el examen.

3) Para averiguar qué técnicas de estudio producen puntuaciones estadísticamente significativas, realiza las siguientes 3 pruebas t por pares:

• Técnica 1 vs. Técnica 2: p=0.0463

• Técnica 1 vs. Técnica 3: p=0.3785

• Técnica 2 vs. Técnica 3: p=0.0114

Quiere controlar la probabilidad de cometer un error de tipo I global (familiar) a \(\alpha\) = 0.05. Decide aplicar una corrección de Bonferroni. ¿Cuál es \(\alpha_{nuevo}\)?

\(\alpha_{nuevo} = \frac{\alpha}{n} = \frac{0.05}{3} = 0.0167\)

4) Procede a comparar p de cada prueba con \(\alpha_{nuevo}\). ¿Cuál es tu observación?

Como el valor p de la técnica 2 vs. técnica 3 es el único valor p<0.0167, concluye que sólo hay una diferencia estadísticamente significativa entre la técnica 2 y la técnica 3.

Advertencia

Bonferroni puede ser demasiado conservador (\(\alpha_{nuevo}\) bastante pequeño), lo que dificulta detectar diferencias que realmente existen (reduciendo su potencia). Pero puede garantizar el error de Tipo I, y es facíl entender la corrección. También lo usamos en comparaciones planeadas.

14.6.2. Prueba Tukey o Tukey HSD

Tukey o Tukey HSD (Honest Significant Difference): mismo tamaño muestral n

Se define el siguiente estadístico:

\[q_N =\frac{X_{i.} - X_{j.}}{\sqrt{\frac{MSE}{N}}}, \qquad df = N = \sum_{i=1}^m n_i-m\]

Tukey-Kramer (también llamado Tukey o Tukey HSD): distintos tamaños muestrales

\[q_N = \frac{X_{i.} - X_{j.}}{\sqrt{MSE \left(\frac{ \frac{1}{n_i}+ \frac{1}{n_j}}{2}\right)}}, \qquad df =N = \sum_{i=1}^m n_i - m\]

donde

\[MSE = \frac{SS_W}{\sum_{i=1}^m n_i - m} = \frac{\sum_{i=1}^m \sum_{j=1}^{n_i} (X_{ij} - X_{i.})^2}{\sum_{i=1}^m n_i - m}\]

En ambos casos \(q_N\) sigue una distribución de rangos student (studendized range distribution) que es similar a t-student (studentized t-test) pero diferente.

Usabmos el \(\alpha\) original (e.g., 0.05). Segimos el proceso de de valor crítico o p-value.

Ejemplo (fuente)

Creamos un conjunto de datos para comparar los rendimientos de 3 modelos. ¿Cuál es \(H_0\)?

#make this example reproducible
set.seed(0)

#create data
df <- data.frame(modelo = rep(c("A", "B", "C"), each = 30),
                 rendimiento = c(runif(30, 0, 3), runif(30, 0, 5), runif(30, 1, 7)))
head(df)
df %>% group_by(modelo) %>% summarise(M=mean(rendimiento), SD=sd(rendimiento))
ggplot(df, aes(x=modelo, y=rendimiento)) + geom_boxplot() + theme(text=element_text(size=14))

A data.frame: 6 x 2

	modelo	rendimiento
	<chr>	<dbl>
1	A	2.6901
2	A	0.7965
3	A	1.1164
4	A	1.7186
5	A	2.7246
6	A	0.6050

A tibble: 3 x 3

modelo	M	SD
<chr>	<dbl>	<dbl>
A	1.584	0.9051
B	2.562	1.2385
C	4.130	1.5683

(Asumimos que los supuestos de ANOVA se mantienen aquí, pero hay que examinar los supuestos en tu análisis de los datos.)

m <- aov(rendimiento~modelo, data=df)
summary(m)

            Df Sum Sq Mean Sq F value  Pr(>F)    
modelo       2   98.9    49.5    30.8 7.5e-11 ***
Residuals   87  139.6     1.6                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

¿Rechazamos \(H_0\)?

Como p<0.05, rechazamos \(H_0\), y decimos que los rendimientos medios de cada modelo no son iguales. Por lo tanto, podemos proceder a realizar la prueba de Tukey para determinar exactamente qué medias de grupo son diferentes.

TukeyHSD(m, conf.level=.95) 

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = rendimiento ~ modelo, data = df)

$modelo
      diff    lwr   upr  p adj
B-A 0.9777 0.1979 1.758 0.0101
C-A 2.5454 1.7656 3.325 0.0000
C-B 1.5677 0.7879 2.347 0.0000

¿Cuál es tu observación?

El valor p indica si existe o no una diferencia estadísticamente significativa entre cada par. Podemos ver que existe una diferencia estadísticamente significativa en los rendimientos entre cada par de modelos al nivel de significación 0.05.

Advertencia

El Tukey es menos conservador que Bonferroni y tiene mejor potencia que Bonferroni, aunque ambos son conservadores (es decir, controlan error de tipo I muy bien) y carecen de potencia estadística.

De las dos, Bonferroni tiene más potencia cuando el número de comparaciones es pequeño, mientras que Tukey es más potente cuando se comprueba un gran número de medias.

En general, Tukey tiene mayor potencia que otras pruebas como Dunn y Scheffé.

Algunas preguntas

Si el test ANOVA es estadísticamente significativo, implica que ¿existe al menos un grupo cuyo comportamiento difiere de los otros?

No necesariamente, podría ocurrir que sea la combinación de grupos que provoque las diferencias (en ANOVA), pero en comparaciones post hoc (e.g., la prueba Tukey o correcdión de Bonferroni) no encontramos diferencias.

Cuando esto ocurre, no hay mucho que podamos hacer, salvo hacer la mejor interpretación que podamos de la inspección de los datos, y ser convenientemente cautos en las conclusiones que saquemos.

Puede resultar tentador realizar un nuevo análisis post hoc utilizando un tipo de prueba diferente. No lo haga. Es una estrategia terrible para hacer estadística, porque este tipo de práctica garantiza el aumento de las posibilidades de encontrar un falso positivo (error de tipo I).

(basado en sugerencia de fuente)

¿Es posible encontrar diferencias siginificativas con las comparaciones múltiples, aunque ANOVA no haya sido significativo?

Si, es posible, debido a que los tests de comparaciones múltiples están mas enfocados (tienen más potencia) que ANOVA.

¿Es útil el resultado del test ANOVA?

Claro, cuando corresponde a la hipótesis en estudio!

A menudo, tus preguntas experimentales están más enfocadas en algunos grupos y sólo quieres comparar ciertos grupos en lugar de todos pares, antes de hacer el experimento/ver los datos. En estos casos, puede hacer comparaciones planeadas. Este fuente tiene más información sobre esto.

14.6.3. Reportar el resultado del ANOVA y comparaciones multiples post hoc

¿Cómo reportamos los resultados? Aquí es un ejemplo:

Según el ANOVA de un factor, hubo un efecto (estadísticamente) significativo del modelo sobre el rendimiento, F(2,87)=30,8, p<0,001. Las comparaciones post hoc mediante la prueba Tukey HSD (o puede sustituirla por la prueba t con corrección de Bonferroni) indicaron que el rendimiento medio del modelo A (M=1,58, SD=0,91) era significativamente diferente del modelo B (M=2,56, SD=1,24; p=0,0101) y del del modelo C (M=4,13, SD=1,57; p<0,001); el rendimiento medio del modelo B era significativamente diferente del del modelo C (p<0,001).

14.7. Variantes de ANOVA y alternativas

Esta clase sólo cubre los fundamento del ANOVA. Para investigaciones o tu proyecto, se puede necesitar modelos más avanzados, y te animo a estudiar los por tu cuenta!

Variantes de ANOVA:

• ANOVA de dos factores (two-way/two-factor ANOVA)

• ANOVA de medidas repetidas (repeated-measure ANOVA)

Alternativas:

• Regresión (veamos esto muy pronto)

• Regresión de efectos mixtos (mixed-effects regression)

14.8. Análisis de asociación en v.a. discretas (o categóricas)

En esta sección se presentan algunos tests que permiten analizar la asociación entre variables aleatorias discretas.

14.8.1. Repaso (independencia)

1) ¿Cómo comprobar si dos eventos \(A\), \(B\) son independientes?

Sean \(A, B \subseteq \Omega\) con \(P(A),P(B) >0\), si cualquiera de las siguientes condiciones se cumple:

1. \(P(A\cap B) = P(A)P(B)\)

2. \(P(A\mid B) = P(A)\)

3. \(P(B\mid A) = P(B)\)

se dice que \(A\) es independiente de \(B\).

2) ¿Cómo comprobar si dos v.a. \(X\), \(Y\) son independientes? (Ejemplo)

Sean \(X\) y \(Y\) dos v.a., discretas o continuas, con distribución de probabilidad conjunta \(f_{X, Y}(x, y)\) y distribuciones marginales \(f_X(x)\) y \(f_Y(y)\), respectivamente. Se dice que las v.a. \(X\) y \(Y\) son estadísticamente independientes si y sólo si

\(f_{X, Y}(x, y) = f_X(x) f_Y(y)\)

para todos los pares \((x,y)\) dentro de sus rangos.

14.8.2. Tabla de contingencia

Consideremos que tenemos observaciones de dos variables discretas o categóricas (o bien variables continuas agrupadas en clases) \(X\) y \(Y\) referidas a una muestra. La tabla de contingencia contiene los números correspondientes a cada pareja de valores (o clases) de ambas variables:

\[\begin{split}\begin{array}{|c||c|c|c|c|c||c|} \hline \bf{X\mid Y}& \bf{d_1} &\cdots &\bf{d_k}& \cdots&\bf{d_s}& \text{total}\\ \hline\hline \bf{c_1} & n_{11} & \cdots & n_{1k}& \cdots& n_{1s} &\bf { n_{1\bullet}}\\ \hline \vdots & \vdots &&\vdots&&\vdots&\vdots\\ \hline \bf{c_h} & n_{h1} & \cdots & n_{hk}& \cdots& n_{hs} & \bf{n_{h\bullet}}\\ \hline \vdots & \vdots &&\vdots&&\vdots&\vdots\\ \hline \bf{c_r} & n_{r1} & \cdots & n_{rk}& \cdots& n_{rs} & \bf{n_{r\bullet}}\\ \hline \hline \text{total} &\bf{n_{\bullet 1}}&\cdots& \bf{n_{\bullet k}}&\cdots&\bf{n_{\bullet s}}&\bf{n}\\ \hline \end{array}\end{split}\]

Cada linea y cada columna corresponden a una submuestra particular. La fila de índice \(h\) es la distribución sobre \( d_1,\ldots,d_s\) de los individuos para los cuales la variable \(X\) toma el valor \(c_h\). La columna de índice \(k\) es la distribución sobre \(c_1,\ldots,c_r\) de los individuos para los cuales las variable \( Y\) toma el valor \(d_k\).

Los totales de renglones y columnas en la tabla \(n_{h\bullet}, n_{\bullet k}\) (etc.) se denominan frecuencias marginales.

Los valores en las celdas internas \(n_{hk}\) (etc.) se denominan frecuencias observadas.

14.8.3. Test de Chi-cuadrado de independencia (para tablas de contingencia)

Tenemos observaciones de dos variables discretas o categóricas (o bien variables continuas agrupadas en clases) \(X\) y \(Y\) referidas a una muestra de distribuciones desconocidas.

Las observaciones de \(X\) son independendes, y las observaciones de \(Y\) son independendes (supuesto de este test).

Nos interesa si \(X\) y \(Y\) son independientes. Tenemos las hipótesis:

\[\begin{split}H_0: \text{X e Y son independientes / no existe relación entre X e Y / X e Y no se afectan mutuamente} \\ H_1: \text{X e Y no son independientes / existe una relación entre X e Y / X e Y se afectan mutuamente}\end{split}\]

Formalmente:

\[\begin{split}H_0: \forall h \in [1,r], k \in [1, s], \qquad p_{hk} = p_h q_k \\ H_1: \exists h \in [1,r], k \in [1, s], \qquad p_{hk} \neq p_h q_k \end{split}\]

donde

• \(p_{hk}\) se refiere a la probabilidad conjunta de que \(X\) sea \(c_h\) y \(Y\) sea \(d_k\)

• \(p_{h}\) se refiere a la probabilidad marginal de que \(X\) sea \(c_h\)

• \(q_{k}\) se refiere a la probabilidad marginal de que \(Y\) sea \(d_k\)

Ya que no conocemos sus distribuciones, estimamos \(p_{hk}, p_{h}, q_{k}\) a partir de la muestra, utilizando la tabla de contingencia:

\[ \hat p_{hk} = \frac{n_{hk}}{n} \qquad \hat p_{h} = \frac{n_{h \bullet}}{n} \qquad \hat q_{k} = \frac{n_{\bullet k}}{n} \]

(Por qué podemos estimar así?)

Asi que bajo \(H_0\), queremos que los valores de \(\hat p_{hk}\) y \(\hat p_{h} \hat q_{k}\) estén cercas:

\[ \frac{n_{hk}}{n} \approx \frac{n_{h \bullet}}{n} \frac{n_{\bullet k}}{n} \]

\[\Leftrightarrow \qquad n_{hk} \text{ (frecuencia observada) } \approx \frac{n_{h \bullet} n_{\bullet k}}{n} \text{ (frecuencia esperada) } \]

\[\Leftrightarrow \qquad O_{hk} \text{ (frecuencia observada) } \approx E_{hk} \text{ (frecuencia esperada) } \]

Definimos un estadístico (v.a.) bajo \(H_0\) que representa la distancia entre ellos para todas combinaciones de valores:

\[TS = \sum\limits_{k=1}^s \sum_{h=1}^r \frac{(O_{hk} - E_{hk})^2}{E_{hk}} = \sum\limits_{k=1}^s \sum_{h=1}^r \frac{\left(n_{hk} - \displaystyle \frac{n_{h \bullet} n_{\bullet k}}{n}\right)^2}{\displaystyle \frac{n_{h \bullet} n_{\bullet k}}{n}} \]

Tenemos un teorema que dice bajo \(H_0\), este estadístico \(TS\) converge hacia una distribución \(\chi^2_{(r-1)(s-1)}\) cuando \( n\) tiende al infinito.

¿Por qué el grado de libertad es \((r-1)(s-1)=rs-r-s+1\)?

Nota

Los grados de libertad se refieren al número de valores que pueden variar libremente en una muestra que se utiliza para estimar las características de una población, dado un número de parámetros estimados antes o un número de restricciones matemáticas.

Usualmente, se encuentran mediante la fórmula \(n-k\), donde \(n\) es el número de observaciones en la muestra y \(k\) es el número de parámetros estimados antes o el número de restricciones.

Es una medida de cuánta información contiene su prueba (test).

Advertencia

Para que la aproximación chi-cuadrado sea válida, cada frecuencia esperada \(E_{hk}\) debe ser de al menos 5 (esta prueba no es válida para muestras pequeñas). Se puede utlizar R para comprobar esto.

Si alguna frecuencia esperada sea menor que 5, consideramos combinar algunas clases o usamos la prueba exacta de Fisher / Fisher’s Exact Test o Monte Carlo simulación (chisq.test(table, simulate.p.value=T)).

Ejemplo (campus y tiempo libre)

Consideremos dos variables categóricas referidas a una muestra de los estudiantes de la UACh, una de las cuales indica el campus, y la otra indica que hace en su tiempo libre. Se trata de analizar si existe una dependencia entre el campus y actividades en su tiempo libre. Supongamos que la tabla de contingencia observada es la siguiente (los números son inventados):

\[\begin{split}\begin{array}{|c||c|c|c|c|} \hline \bf{X\mid Y}& \text{Deporte} & \text{Estudio} &\text{Trabajo} & \text{Total} \\ \hline\hline \text{Campus Isla Teja} & 68 & 56 & 32 & 156\\ \hline \text{Campus Miraflores} & 52 & 72 & 20 & 144\\ \hline \text{Total} &120 & 128 & 52 & 300\\ \hline \end{array}\end{split}\]

(Asumimos que cada frecuencia esperada \(E_{hk} \geq 5\), pero tienes que examinar esto en tu análisis de los datos.)

¿Cuál es \(H_0\) aquí?

En tu tarea, tienes que desarrollar el test. Aquí, mostramos el resultado de R:

#usando el test predifinido en R
table <- rbind(c(68, 56, 32), c(52, 72, 20))
chisq.test(table, correct=FALSE)

	Pearson's Chi-squared test

data:  table
X-squared = 6.4, df = 2, p-value = 0.04

¿Se recharza \(H_0\) o no? ¿Qué significa este resultado?

Se recharza \(H_0\) en favor de dependencia entre campus y actividades en tiempo libre. Es decir, tenemos evidencia para apoyar que el campus afecta las actividades en tiempo libre.

¿Cómo reportar esta prueba formalmente?

Se realizó una prueba chi-cuadrado de independencia para examinar la relación entre [variable 1] y [variable 2]. Hubo [o no hubo] una relación significativa entre las dos variables, \(\chi^2\)(df, N) = [valor \(\chi^2\)], p = [valor p].

(hay flexibilidad en la presentación de informes que varía un poco).

14.8.4. Test de McNemar-Bowker (o Bowker)

Este test se utiliza para analizar la simetría en una tabla de contingencia construida a partir de dos muestras dependientes \(X, Y\) (donde las observaciones de cada muestra son independentes), típicamente pre y post un tratamiento.

Las hípotesis de este test es:

\[\begin{split}H_0: \forall i, j \in [1,k], i \neq j, \qquad p_{ij} = p_{ji} \\ H_1: \exists i, j \in [1,k], i \neq j, \qquad p_{ij} \neq p_{ji}\end{split}\]

donde

• \(p_{ij}\) se refiere a la probabilidad de que \(X\) sea \(i\) y \(Y\) sea \(j\) (usamos \(i\) o \(j\) para referir a un valor de una clase o categoría, aunque strictamente debemos escribir \(c_i\) o \(c_j\))

• \(p_{ji}\) se refiere a la probabilidad de que \(X\) sea \(j\) y \(Y\) sea \(i\)

El estadístico del test de Bowker se construye como:

\[ B = \sum_{j>i} \frac{(n_{ij} - n_{ji})^2}{n_{ij} + n_{ji}}\]

De manera que bajo la hipótesis nula de simetría en la tabla de contingencia, es decir que no hay efecto del tratamiento, se cumple asintóticamente que:

\[B \sim \chi^2_{\frac{k(k-1)}{2}}\]

Si se rechaza la hipótesis nula (cuándo?), y además se observan más efectivos sobre la diagonal (upper off-diagonal) que bajo la diagonal (lower off-diagonal) se puede concluir que hay un efecto positivo del tratamiento en el criterio en análisis.

Advertencia

• Para que la aproximación chi-cuadrado sea válida, se recomienda que \(n_{ij} + n_{ji} \geq 10, \forall i, j, i \neq j\) (o \( n_{ij} + n_{ji} \geq 25\) según otro fuente). Si no se cumple, consideramos combinar algunas clases o usamos el mid-P McNemar (binomial) test.

• Si \(n_{ij} + n_{ji} = 0\) en algún par, lo ignoramos en la calculación de B.

Ejemplo

Considere los datos de 170 estudiantes sometidos a una evaluación sobre sus habilidad de cohesión en la redacción de textos, antes y después de una intervención didáctica. La evaluación utiliza una escala de 4 niveles desde no logrado a totalmente logrado. Se pide evaluar si la intervención ha tenido efecto o no en la habilidad estudiada.

(Asumimos que \(n_{ij} + n_{ji} \geq 10, \forall i, j, i \neq j\) (o \( n_{ij} + n_{ji} \geq 25\) según otro fuente), pero tienes que examinar esto en tu análisis de los datos.)

#fila (i): pre-evaluación
#columna (j): post-evaluación
tabla <- rbind(c(0,6,9,6), c(6,6,9,21), c(9,12,18,55), c(4,0,3,6))
print(tabla)
B = 0
efsd = efbd = 0
for (i in 1:3){
    for (j in (i+1):4){
        B = B + (tabla[i,j] - tabla[j,i])^2 / (tabla[i,j] + tabla[j,i])
        efsd = efsd + tabla[i,j]
        efbd = efbd + tabla[j,i]
    }
}
cat("B:", B, "\n")
cat("sobre diag vs. bajo diag:", c(efsd, efbd), "\n")

k = 4*3/2
q = qchisq(0.95, k)
print(q)

p <- 1-pchisq(B, k)
print(p)

     [,1] [,2] [,3] [,4]
[1,]    0    6    9    6
[2,]    6    6    9   21
[3,]    9   12   18   55
[4,]    4    0    3    6
B: 68.44926 
sobre diag vs. bajo diag: 106 34 
[1] 12.59159
[1] 8.500978e-13

Podemos comprobar los resultados anteriores con mcnemar.test en R:

##usando el test predefinido en R
mcnemar.test(tabla)

	McNemar's Chi-squared test

data:  tabla
McNemar's chi-squared = 68.449, df = 6, p-value = 8.501e-13

¿Se rechaza \(H_0\)? ¿Qué significa este resultado?

Se rechaza \(H_0\) de simetría en favor de un efecto positivo de la intervención didáctica.